home *** CD-ROM | disk | FTP | other *** search
/ Collection of Tools & Utilities / Collection of Tools and Utilities.iso / c / fftsing.zip / SING.C < prev    next >
C/C++ Source or Header  |  1992-04-15  |  18KB  |  844 lines
  1. /**************************************************************************
  2.  
  3.     Javier Soley, Ph. D,   FJSOLEY @UCRVM2.BITNET
  4.     Escuela de Física y Centro de Investigaciones Geofísicas
  5.     Universidad de Costa Rica
  6.  
  7. ***************************************************************************/
  8.  
  9. /*    Computes the DISCRETE FOURIER TRANSFORM of very long data series.
  10.  *   Restriction: the data has to fit in conventional memory.
  11.  *
  12.  *   Compile in compact or large models ===> Pointers must be FAR
  13.  *   Two huge pointers are used to access the real and imaginary
  14.  *   parts of the transform without 64k wrap around. 
  15.  *
  16.  *   This functions are translations from the fortran program in
  17.  *
  18.  *   R. C. Singleton, An algorithm for computing the mixed radix fast
  19.  *   Fourier transform 
  20.  *
  21.  *    IEEE Trans. Audio Electroacoust., vol. AU-17, pp. 93-10, June 1969.
  22.  *   Some features are:
  23.  *
  24.  *        1-) Accepts an order of transform that can be factored not only
  25.  *            in prime factors such 2 and 4, but also including odd factors
  26.  *            as 3, 5, 7, 11, etc.
  27.  *        2-) Generates sines and cosines recursively and includes
  28.  *            corrections for truncation errors.
  29.  *        3-) The original subroutine accepts multivariate data. This 
  30.  *            translation does not implement that option (because I
  31.  *            do not needed right now).
  32.  *
  33.  *    Singleton wrote his subroutine in Fortran and in such a way that it 
  34.  *    could be ported allmost directly to assembly language. I transcribed
  35.  *   it to C with little effort to make it structured. So I apologize to
  36.  *   all those C purists out there!!!!!!!!
  37.  *
  38.  */
  39.  
  40.                 /*  Version 2.0 March/30/92 */
  41. /* Includes */
  42.  
  43. #include <stdlib.h>
  44. #include <stdio.h>
  45. #include <math.h>
  46.  
  47. /* Defines */
  48.  
  49. #define    TWO_PI    ((double)2.0 * M_PI)
  50. #define   MAXF  23
  51. #define   MAXP  209
  52.  
  53. #define   MY_SIZE     double   /* or float or long double */    
  54. #define   MY_TYPE     huge     /* or far if data fits in two 64 K segments */
  55.  
  56.  
  57. /* Globals */
  58.  
  59. long nn, m, flag, jf, jc, kspan, ks, kt, nt, kk, i;
  60. MY_SIZE   c72, s72, s120, cd, sd, rad, radf, at[23], bt[23];
  61. long nfac[23];   
  62. int inc;
  63. long np[MAXP];
  64.  
  65.  
  66. /* The functions */
  67.  
  68. void radix_2( a, b)
  69.  
  70. MY_SIZE MY_TYPE *a;
  71. MY_SIZE MY_TYPE *b;   
  72.  
  73. { long k1, k2;
  74.   MY_SIZE ak, bk, c1, s1;
  75.     kspan >>= 1;
  76.     k1 = kspan +2;
  77.  
  78.     do {
  79.         do  {
  80.           k2 = kk + kspan;
  81.           ak = a[k2-1];
  82.           bk = b[k2-1];
  83.           a[k2-1] = a[kk-1] -ak;
  84.           b[k2-1] = b[kk-1] -bk;
  85.           a[kk-1] += ak;
  86.           b[kk-1] += bk;
  87.           kk = k2 + kspan;
  88.         } while ( kk <= nn);
  89.        kk = kk - nn;
  90.     } while ( kk <= jc);
  91.  
  92.  
  93.     if ( kk > kspan) flag = 1; 
  94.     else
  95.     {
  96.     do {
  97.         c1 = 1.0 - cd;
  98.         s1 = sd;
  99.         do {
  100.             do  {
  101.                 do {
  102.                     k2 = kk + kspan;
  103.                     ak = a[kk-1]- a[k2-1];
  104.                     bk = b[kk-1]- b[k2-1];
  105.                     a[kk-1] += a[k2-1];
  106.                     b[kk-1] += b[k2-1];
  107.                     a[k2-1]  = c1*ak - s1*bk;
  108.                     b[k2-1]  = s1*ak + c1*bk;
  109.                     kk = k2 + kspan;
  110.                 } while ( kk < nt );
  111.             k2 = kk - nt;
  112.             c1 = -c1;
  113.             kk = k1 - k2;
  114.             } while ( kk > k2 );
  115.         ak = c1- (cd*c1+sd*s1);
  116.         s1 = (sd*c1-cd*s1) +s1;
  117.  
  118.      /***** Compensate for truncation errors   *****/
  119.  
  120.         c1 = 0.5/(ak*ak+s1*s1)+0.5;
  121.         s1 *= c1;
  122.         c1 *= ak;
  123.         kk += jc;
  124.         } while ( kk < k2);
  125.         k1 = k1 + inc + inc;
  126.         kk = (k1- kspan) /2 + jc;
  127.        } while ( kk <= jc + jc );
  128.     }
  129. }
  130.  
  131. void radix_4( isn, a, b)
  132.  
  133. MY_SIZE MY_TYPE *a;
  134. MY_SIZE MY_TYPE *b;
  135. int  isn;
  136. {
  137.     long    k1, k2, k3;
  138.     MY_SIZE  akp, akm, ajm, ajp, bkm, bkp, bjm, bjp;
  139.     MY_SIZE  c1, s1, c2, s2, c3, s3;
  140.     kspan /= 4;
  141. cuatro_1: 
  142.     c1 = 1.0;
  143.     s1 = 0;
  144.     do {
  145.         do {
  146.             do {
  147.                 k1  =  kk + kspan;
  148.                 k2  =  k1 + kspan;
  149.                 k3  =  k2 + kspan;
  150.                 akp =  a[kk-1] + a[k2-1];
  151.                 akm =  a[kk-1] - a[k2-1];
  152.                 ajp =  a[ k1-1] + a[k3-1];
  153.                 ajm =  a[ k1-1] - a[k3-1];
  154.                 a[kk-1] = akp + ajp;
  155.                 ajp = akp - ajp;
  156.                 bkp = b[kk-1] + b[k2-1];
  157.                 bkm = b[kk-1] - b[k2-1];
  158.                 bjp = b[k1-1] + b[k3-1];
  159.                 bjm = b[k1-1] - b[k3-1];
  160.                 b[kk-1] = bkp + bjp;
  161.                 bjp = bkp - bjp;
  162.                 if ( isn < 0) goto cuatro_5;
  163.                 akp = akm - bjm;
  164.                 akm = akm + bjm;
  165.                 bkp = bkm + ajm;
  166.                 bkm = bkm - ajm;
  167.                 if (s1 == 0.0) goto cuatro_6;
  168.     cuatro_3:        a[ k1-1] = akp*c1 - bkp*s1;
  169.                 b[ k1-1] = akp*s1 + bkp*c1;
  170.                 a[ k2-1] = ajp*c2 - bjp*s2;
  171.                 b[ k2-1] = ajp*s2 + bjp*c2;
  172.                 a[ k3-1] = akm*c3 - bkm*s3;
  173.                 b[ k3-1] = akm*s3 + bkm*c3;
  174.                 kk = k3 + kspan;
  175.                 }  while ( kk <= nt);
  176.                 
  177.     cuatro_4: 
  178.         c2 = c1 - (cd*c1 + sd*s1);
  179.         s1 = (sd*c1 - cd*s1) + s1;
  180.  
  181.     /***** Compensate for truncation errors *****/
  182.  
  183.         c1 = 0.5 / (c2*c2 + s1*s1) +0.5;
  184.         s1 = c1 * s1;
  185.         c1 = c1 * c2;
  186.         c2 = c1*c1 - s1*s1;
  187.         s2 = 2.0 * c1 *s1;
  188.         c3 = c2*c1 - s2*s1;
  189.         s3 = c2*s1 + s2*c1;
  190.         kk = kk -nt + jc;
  191.         } while ( kk <= kspan);
  192.         
  193.     kk = kk - kspan + inc;
  194.     if ( kk <= jc)  goto cuatro_1;
  195.     if ( kspan == jc)  flag =1;
  196.     goto out;
  197. cuatro_5:
  198.     akp = akm + bjm;
  199.     akm = akm - bjm;
  200.     bkp = bkm - ajm;
  201.     bkm = bkm + ajm;
  202.     if (s1 != 0.0) goto cuatro_3;
  203. cuatro_6:
  204.     a[k1-1] = akp;
  205.     b[k1-1] = bkp;
  206.     b[k2-1] = bjp;
  207.     a[k2-1] = ajp;
  208.     a[k3-1] = akm;
  209.     b[k3-1] = bkm;
  210.     kk = k3 + kspan;
  211.     } while ( kk <= nt);
  212.     goto cuatro_4;
  213. out: 
  214.     s1 = s1 + 0.0;
  215. }
  216.  
  217.  
  218.  
  219.  
  220.     /* Find prime factors of n */
  221.  
  222. void fac_des( n )
  223. long n;
  224. {
  225.     long k, j, jj, l;
  226.     k  = n;
  227.     m = 0;
  228.     while ( k-(k / 16)*16 == 0 ) {
  229.         m++;
  230.         nfac[m-1] = 4;
  231.         k /= 16;
  232.     } 
  233.     j  = 3;
  234.     jj = 9;
  235.     do {
  236.          while (k % jj == 0) {
  237.              m++;
  238.              nfac[m-1] = j;
  239.              k /= jj;
  240.          }
  241.     j += 2;
  242.     jj = j * j;
  243.     } while ( jj <= k);
  244.     if (k <= 4) {
  245.         kt = m;
  246.         nfac[m] = k;
  247.         if (k != 1) m++;
  248.     }
  249.     else {
  250.         if (k-(k / 4)*4 == 0) {
  251.             m++;
  252.             nfac[m-1] = 2;
  253.             k /= 4;
  254.         }
  255.         kt = m;
  256.         j = 2;
  257.         do {
  258.             if (k % j == 0 ) {
  259.                 m++;
  260.                 nfac[m-1] = j;
  261.                 k /= j;
  262.             }
  263.         j = ((j+1)/ 2)*2 + 1;
  264.         } while ( j <= k);
  265.     }
  266.     if (kt != 0) {
  267.         j = kt;
  268.         do {
  269.            m++;
  270.            nfac[m-1] = nfac[j-1];
  271.            j--;
  272.         } while ( j != 0);
  273.     }
  274. }    
  275.  
  276.     /* Permute the results to normal order  */
  277.  
  278.  
  279. void permute( ntot, n, a, b)
  280.  
  281. long ntot, n;
  282. MY_SIZE MY_TYPE *a;
  283. MY_SIZE MY_TYPE *b;
  284.  
  285.  
  286.  
  287. {
  288.     long  k, j, k1, k2, k3, kspnn, maxf;
  289.     
  290.     MY_SIZE ak, bk;
  291.     long  ii, jj;
  292.     maxf = MAXF;
  293.     np[0] = ks;
  294.     if (kt != 0) {
  295.         k = kt +kt +1;
  296.         if (m < k) k--;
  297.         j = 1;
  298.         np[k] = jc;
  299.         do {
  300.             np[j]   = np[j-1] / nfac[j-1];
  301.             np[k-1] = np[k]   * nfac[j-1];
  302.             j++;
  303.             k--;
  304.         } while (j < k);
  305.         k3 = np[k];   
  306.         kspan = np[1];
  307.         kk = jc+1;
  308.         k2 = kspan + 1;  
  309.         j = 1;
  310.  
  311.  /***** Permutation of one dimensional transform *****/
  312.  
  313.         if (n == ntot) {
  314.             do {
  315.                 do {
  316.                     ak      = a[kk-1];
  317.                     a[kk-1] = a[k2-1];
  318.                     a[k2-1] = ak;
  319.                     bk      = b[kk-1];
  320.                     b[kk-1] = b[k2-1];
  321.                     b[k2-1] = bk;
  322.                     kk     += inc;
  323.                     k2     += kspan;
  324.                 } while ( k2 < ks);
  325. ocho_30:          do {
  326.                     k2 -= np[j-1];
  327.                     j++;
  328.                     k2 += np[j];
  329.                 } while (k2 > np[j-1]);
  330.                 j = 1;
  331. ocho_40:          j = j + 0;
  332.             } while (kk < k2);
  333.             kk += inc;
  334.             k2 += kspan;
  335.             if (k2 < ks)  goto ocho_40;
  336.             if (kk < ks)  goto ocho_30;
  337.             jc = k3;
  338.         }
  339.         else {        /* Permutation for multiple transform  */
  340. ocho_50:
  341.             do {
  342.                 do {
  343.                     k = kk + jc;
  344.                     do {
  345.                         ak = a[kk-1];    
  346.                         a[kk-1] = a[k2-1];
  347.                         a[k2-1] = ak;                            
  348.                         bk = b[kk-1];
  349.                         b[kk-1] = b[k2-1];
  350.                         b[k2-1] = bk;
  351.                         kk += inc;
  352.                         k2 += inc;
  353.                     } while ( kk < k);
  354.                     kk = kk + ks - jc;
  355.                     k2 = k2 + ks - jc;
  356.                 } while (kk < nt);
  357.                 k2 = k2 - nt + kspan;
  358.                 kk = kk - nt + jc;
  359.             } while (k2 < ks);
  360.             do {
  361.                 do {
  362.                     k2 -= np[j-1];
  363.                     j++;
  364.                     k2 += np[j];
  365.                 } while (k2 > np[j-1]);
  366.                 j =1;
  367.                 do {
  368.                     if ( kk < k2 ) goto ocho_50;
  369.                     kk +=jc;
  370.                     k2 += kspan;
  371.                 } while (k2 < ks);
  372.             } while (kk < ks);
  373.             jc = k3;
  374.         }
  375.     }
  376.  
  377.     if ( (2*kt +1) < m) {
  378.         kspnn = np[kt];
  379.                 /* Permutation of square-free factors of n */
  380.         j = m - kt;
  381.         nfac[j] = 1;
  382.         do {
  383.             nfac[j-1] *= nfac[j];
  384.             j--;
  385.         } while (j != kt);
  386.         kt++;
  387.         nn = nfac[kt-1] -1;
  388.         if (nn > MAXP) {
  389.             printf("product of square free factors exceeds allowed limit\n");
  390.             exit(2);
  391.         }
  392.         jj =0;
  393.         j=0;
  394.         goto nueve_06;
  395. nueve_02:
  396.         jj -= k2;
  397.         k2 = kk;
  398.         k++;
  399.